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GENERALIZED  HYDRODYNAMICS  WITH  VISCOPLASTICITY 
FOR  CHANNELING  IN  SATURATED  SAND 


INTRODUCTION 

An  important  issue  in  littoral  warfare  is  the  clearance  of  obstacles  and  mines  to  enable  a  beach 
landing.  One  possible  solution  to  this  problem  is  to  use  explosives  to  form  craters  and  channels  in  the 
sandy  soil  [1,  2].  The  simulation  of  this  phenomenon  offers  a  great  challenge  to  computational  models, 
not  only  because  of  violent  ffee-surface  motions  but  also  due  to  the  complex  interaction  of  the  water  with 
the  sand  bottom.  This  report  describes  an  attempt  to  model  this  problem,  using  a  viscoplastic  description 
for  saturated  sand,  within  a  previously  developed  method  for  simulating  violent  ffee-surface  motions. 

In  a  single-charge  underwater  explosion,  a  gas  bubble  forms  immediately  after  the  detonation  shock 
wave.  This  bubble  is  initially  composed  of  high-pressure  gas  that  causes  the  bubble  to  grow 
spherically  [3].  As  the  bubble  radius  increases,  its  pressure  drops  (adiabatically).  When  the  bubble 
pressure  drops  below  the  ambient  pressure,  its  radius  will  begin  to  decelerate  and  the  bubble  will 
eventually  attain  a  maximum  volume.  At  this  time,  the  bubble  pressure  is  at  a  minimum  and  the  bubble 
will  begin  to  collapse.  Except  at  very  great  depths  or  in  cases  when  any  boundaries  are  sufficiently  far 
from  the  charge  and  gravity  is  absent,  the  bubble  will  usually  not  remain  spherical  during  its  collapse. 
Indeed,  at  relatively  shallow  depths  or  cases  when  a  boundary  is  near  the  charge,  the  bubble  will  lose  its 
spherical  shape  before  it  attains  its  maximum  size.  Bubble  jetting  usually  occurs  during  the  collapse 
phase.  This  can  be  predicted  either  by  theory  [4]  or  by  direct  simulation.  For  example,  in  a  shallow  depth 
explosion  (which  does  not  vent),  a  plume  of  water  is  ejected  upward  above  the  bubble  as  the  bubble 
attains  its  maximum  volume.  At  the  same  time,  another  jet  of  water  flows  downward  from  the  plume 
through  the  bubble  as  it  begins  its  collapse  [4,  5,  6]. 

Consider  the  case  when  the  charge  is  placed  in  shallow  water  above  a  sandy  bottom.  A  crater  can  be 
expected  to  form  as  the  high-pressure  gas  causes  the  bubble  to  expand,  pushing  the  water  upward  and 
saturated  sand  downward.  If  the  bubble  does  not  vent,  the  downward  moving  jet  will  impact  the  bottom  of 
the  sand  crater  and  cause  further  excavation.  At  later  times  the  sediment  from  the  sand-water  mixture 
refills  part  of  the  crater.  The  dynamics  of  the  sand-water  slurry  during  the  violent  collapse  of  the  bubble 
and  subsequent  re-expansion  are  not  well  understood. 

As  can  be  expected,  the  dynamics  and  the  final  crater  size  will  depend  on  the  properties  of  the 
saturated  sand.  This  dependence  was  noted  in  Ref.  7,  Ch.  6  in  which  empirical  laws  were  given  for  crater 
radius  as  a  function  of  the  charge  depth  (height  of  water).  Due  to  the  variance  in  the  types  of  bottoms,  the 
empirical  laws  for  crater  radius  were  stated  to  have  a  general  prediction  accuracy  of  +30%.  It  was  also 
reported  in  Ref.  7  that  the  crater  depths  have  an  even  greater  variance  with  soil  type.  The  most  dramatic 
example  of  this  dependence  is  the  experimental  results  using  similar  charge  configurations  at  different 
sites  as  indicated  in  Refs.  8  and  9.  For  example,  explosions  in  shallow  water  conducted  at  Port  Wakefield, 
Australia  [8],  always  produced  craters  or  channels  in  the  sand.  However,  shallow  water  explosions  over 
the  sand  in  the  Eglin  Air  Force  Base  Test  Pond  [9]  caused  the  sand  to  be  raised,  forming  a  berm. 
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In  previous  work  [5,  6,  10-13],  the  authors  and  colleagues  have  developed  a  computational  capability 
for  free  surface  flows  in  an  incompressible  fluid  based  on  a  generalized  hydrodynamics  model.  The  code 
BUB2D  has  been  used  to  simulate  two-dimensional  (2-D)  or  axially  symmetric  problems,  and  BUB3D 
simulates  fully  three-dimensional  (3-D)  flows.  The  generalized  model  provides  a  method  for  density  and 
momentum  redistribution,  in  such  a  way  that  the  density  constraint  on  the  incompressible  liquid  is 
satisfied,  and  energy  is  nonincreasing.  The  energy  property  of  this  method  makes  it  unique  among  volume 
of  fluid  (VOF)  type  procedures,  which  usually  simply  truncate  any  excesses  in  density  without  any 
guarantee  that  mass  is  conserved  or  energy  is  nonincreasing  [14,  15].  In  addition  to  the  improved  stability 
provided  by  the  unique  redistribution  method,  stable  and  accurate  procedures  are  employed  for  the  other 
steps  in  the  solution  process.  For  example,  a  second-order  Godunov  method  is  used  for  the  treatment  of 
the  nonlinear  convective  terms,  and  an  approximate  pressure  projection  method  is  employed  for  the 
determination  of  the  pressure,  so  that  the  flow  is  divergence  free  in  the  liquid  region. 

While  these  codes  have  been  applied  to  many  flow  situations,  the  most  relevant  to  surf-zone 
explosive  cratering  problems  are  the  examples  of  shallow  depth  explosions  in  a  single  incompressible 
species  (water)  [5,  6].  The  presence  of  saturated  sand  complicates  the  situation,  primarily  due  to  its  dual 
dynamical  properties.  For  example,  if  a  bed  of  sand  is  shaken  (or  driven  with  enough  force),  it  becomes 
fluidized  [16-18].  However,  sand  (particularly,  wet  sand)  can  obviously  behave  as  a  solid,  a  fact  that 
anyone  who  has  built  a  sandcastle  can  attest  to.  This  dual  property  suggests  that  some  form  of  “plasticity” 
is  required  for  an  appropriate  model. 

Traditional  continuum  models  for  granular  flows  are  usually  based  on  the  Mohr-Coulomb  criterion 
for  failure.  This  criterion  states  that  yielding  will  occur  on  a  plane  element  when 

151  =  c  +  T  tan(0) , 

where  S  is  the  shear  stress,  c  is  the  cohesion,  T  is  the  normal  stress,  and  <p  is  the  internal  angle  of 
friction  for  the  granular  material.  Constitutive  laws,  which  have  incorporated  this  criterion,  have  been 
presented  in,  for  example,  Refs.  19  to  21.  Alternately,  a  Drucker-Prager  discontinuity  condition  was  used 
in  place  of  the  Mohr-Coulomb  criterion  in  Refs.  22  and  23.  This  condition  was  claimed  to  have 
computational  advantages  over  the  Mohr-Coulomb  criterion,  while  producing  similar  results  [22]. 

For  the  present  study,  however,  we  have  chosen  a  Bingham  plastic  model  for  the  saturated  sand.  This 
model  provides  a  simple  viscoplastic  description,  consisting  of  only  a  yield  stress  and  a  viscosity  during 
yielding  [24].  It  is  typically  used  to  describe  the  rheology  in  fluids  with  a  large  amount  of  suspended 
solids  [25]  and  has  been  used  in  the  description  of  debris  flow  [26].  Furthermore,  for  channel-flow 
problems,  models  using  the  Mohr-Coulomb  criterion  have  been  shown  to  exhibit  a  plug  region  in  the 
center  of  the  channel,  similar  to  that  of  a  Bingham  fluid  [19].  More  realistic  descriptions  of  rate  dependent 
stress  models  can  be  easily  incorporated  into  our  numerical  algorithm,  but  this  report  only  discusses  the 
capabilities  and  limitations  of  applying  the  Bingham  model. 

This  report  begins  by  providing  a  brief  outline  of  the  generalized  hydrodynamics  model  and  its 
extension  (in  particular,  the  redistribution  phase)  to  multiple  species.  A  description  of  material  models, 
including  plasticity,  is  given  in  the  next  section.  Particular  emphasis  is  placed  on  the  use  of  non- 
Newtonian  viscosity,  in  particular,  a  Bingham  plastic  model  for  the  saturated  sand.  The  details  of 
implementing  an  approximation  to  the  Bingham  plastic  model  within  the  context  of  generalized 
hydrodynamics  are  also  presented.  Validations  are  shown  for  benchmark  problems,  including 
comparisons  to  an  analytic  solution  of  flow  in  a  circular  cylinder.  Finally,  a  comparison  is  made  from 
videos  taken  from  an  experiment  and  computations  using  various  values  for  the  yield  stress  within  the 
Bingham  model  for  sand. 
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MULTIPLE  SPECIES  FORMULATIONS 

In  this  section,  a  description  of  our  generalized  hydrodynamics  model  is  presented  together  with  the 
modifications  required  for  the  treatment  of  multiple  incompressible  species.  We  start  with  the  simple 
inviscid  case;  the  formulation  and  implementation  of  viscosity  are  described  in  detail  in  the  succeeding 
sections. 

Two  Species  Model 

For  simplicity,  first  consider  the  two-species  case  where,  for  example,  the  species  are  saturated  sand 
and  water.  For  the  remainder  of  this  paper  we  will  refer  to  “saturated  sand”  as  simply  “sand.”  Let  p  be  the 
mixture  density, 


p  =  pA+pJw *  W 

where  ps  is  the  density  of  sand,  pw  is  the  density  of  water,  and  </>s  >  0,  (j)w>  0  are  the  volume  fractions  of 
sand  and  water,  respectively.  For  an  inviscid  fluid  mixture  in  which  the  two  species  move  with  a  common 
velocity  (without  any  material  strength),  the  Conservation  of  Momentum  equation  is 

(pu),  +  V  •  (puu)  =  -VP  -  pgk,  (2) 

where  u  =  (u,v,wj  is  the  fluid  (mixture)  velocity,  P  is  the  pressure,  g  is  the  acceleration  due  to  gravity, 

and  k  is  the  unit  vector  in  the  direction  opposite  to  the  gravity.  (In  the  actual  formulation  of  the 
Generalized  Formulation  of  Hydrodynamics  considered  here,  the  pressure  gradient  does  not  explicitly 
appear  in  the  momentum  conservation.  Instead,  it  is  added  in  implicitly  as  a  Lagrange  multiplier  used  to 
satisfy  constraints  on  the  density  and  the  effect  of  these  constraints  on  the  conservation  of  mass  equation. 
For  a  more  detailed  discussion  on  this,  see  Refs.  10  to  13.)  The  Conservation  of  Mass  may  be  expressed 
as  conservation  of  each  species  volume;  namely, 


(<Ps)t  +V*(0su)  =  O,  (3) 

(<t>w),  +V*(0wu)  =  O.  (4) 

Equations  (3)  and  (4)  are  solved  with  the  constraints: 

(5) 

^  1,  (6) 

+  !•  (7) 


Note,  the  only  actual  constraint  is  Eq.  (7)  since  Eq.  (7)  together  with  the  non-negativity  of  the  volume 
fractions  yields  Eqs.  (5)  and  (6).  Also,  the  non-negativity  of  these  variables  is  preserved  by  the  solution  of 
Eqs.  (3)  and  (4),  so  that  <j)>  0  need  not  be  explicitly  imposed  as  a  constraint.  Note,  that  Eqs.  (3)  and  (4) 

yield 


V*u  =  0  when  (j)w  +  (j)s  =1, 


(8) 


which  is  the  usual  divergence  free  condition  for  incompressible  flow.  Also  Eqs.  (1),  (3),  and  (4)  yield 
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P,  +V*(pu)  =  0,  (9) 

which  is  the  usual  mass  conservation  equation. 

This  constrained  system  of  conservation  laws  can  be  approximated  using  the  Generalized 
Formulation  of  Hydrodynamics  for  Multiple  Species.  For  single  species,  this  formulation  was  first 
proposed  in  Ref.  10.  This  formulation  may  be  stated  as  an  algorithm  that  advances  the  flow  at 

time  f  to  the  flow  (pn+1 ,  u”+l )  at  time  t"+1  =  t”  +  T . 

The  first  step  is  to  approximate  the  solution  to  Eqs.  (2)  to  (4)  without  the  constraints  Eqs.  (5)  to  (8) 
and  without  the  pressure  term  in  the  momentum  equations.  This  solution  can  be  obtained  by 
approximating  the  coupled  Eqs.  (2)  and  (9)  (which  are  the  same  as  for  a  single  incompressible  species) 
and  including  either  Eqs.  (3)  or  (4),  (say  Eq.  (3)).  We  denote  this  flow  as  (p,u)  and  let  <j)s  denote  the 
intermediate  “sand”  volume  fraction.  Then 
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and 


h.=H"  if  ^>0 

1  0,  if  Hs  =0 

h-h-tf,  if  Hw> 0 
"  1  0,  if  Hw=  0’ 


(14) 

(15) 


with  the  redistributed  volume  fractions  given  by 

<t>f=l+V2Hs  (16) 

and 

o7) 

Applying  the  maximum  principle  to  Eq.  (14)  yields  V2//s  <  1  —  .  This,  together  with  Eq.  (16),  implies 

that  <j)fl  <1.  Similarly,  Eqs.  (15)  and  (17)  ensure  that  <pf'  +  <pfl  <1.  Finally,  Eqs.  (1),  (13),  (16),  and 
(17)  yield 

p,+,=ftr+pi:+1=p+v2^  as) 


which  is  consistent  with  Eq.  (10). 

The  final  step  carries  the  solution  from  (p"+1,u)  to  (pn+',u"+1)  so  that  the  constraint  Eq.  (8)  is 

satisfied.  This  is  done  by  projecting  the  velocity  onto  its  divergence  free  subspace  in  the  domain  where 
<Pw  +(j)s  =  1.  This  projection  involves  the  addition  of  the  gradient  of  the  pressure  (used  here  as  a  Lagrange 
multiplier)  to  the  momentum.  That  is 

p"+V+1=pn+1u-zVP,  (19) 

where  P  solves  the  equation 

*V*-L-VP  =  V«u  in  Dn+l,  (20) 

P 

where  Dn+l  is  the  “incompressible”  region  where  (j)w  +  (f>s  =  1 .  In  practice,  this  region  is  defined  to  be  the 
region  where  </>w  +  <j)s  >  1  -  £  for  some  prescribed  value  of  e.  Outside  of  this  region,  the  pressure  is 
assumed  to  be  uniform  within  each  component.  For  example,  in  the  “air”  region,  the  pressure  is  set  to  the 
ambient  pressure  PA-  Inside  a  bubble,  the  pressure  is  determined  using  an  adiabatic  law.  Note,  that  Eqs. 
(19)  and  (20)  yield 


V  •  u”+1 


=  V  •  u  -  rV  • 


VP  =  0 


in  D" 


(21) 


General  m-Species  Case 


In  the  case  when  there  are  m  species  and 
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p  =  £pA. 

1=1 


(22) 


the  transport  species  equations  are 

(0,),  +  V •  (0,-u)  =  0 ,  for /  =  1,  ...,/m. 


(23) 


771 

The  successive  constraint  procedure  for  imposing  0,  <  1  is  accomplished  by  setting 

7=1 


where  ft  solves 


Then 


V2ft  =  J 


i-1 


v  Ht>0  f  •  , 

"  ;  ,  for  i  =  1, m. 

0,  if  H(=  0 


P"+i  =  £p,0r'  =  P  +  £p,V2tf,  =  p  +  V2H, 


7=1 


(24) 


(25) 


(26) 


where 


H=±PiHt.  (27) 

7=1 

As  with  the  two-species  formulation,  Eq.  (11)  is  used  for  the  momentum  redistribution  with  H  given  by 
Eq.  (27),  and  energy  cannot  increase.  The  pressure  projection  is  performed  exactly  as  before,  following 
Eqs.  (19)  to  (21). 

VISCOSITY  AND  PLASTIC  DEFORMATION 

In  this  section,  we  discuss  some  material  models  including  both  Newtonian  and  non-Newtonian 
viscosity,  and  elastic-plastic  deformation.  The  implementation  of  an  “exponential”  non-Newtonian 
viscous  model,  which  can  closely  approximate  a  Bingham  plastic  model,  is  discussed  in  some  detail. 

Material  Models 

Consider  the  momentum  equation  with  the  inclusion  of  viscous  and/or  material  stresses, 

(pu), +V*(puu)  =  -V»o-pgk,  (28) 


with 


a  =  PI-T, 


(29) 
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where  P  is  the  pressure,  I  is  the  3  X  3  identity  matrix  (tensor),  and  T  is  the  stress  tensor.  In  the  case 
whenT  =  T(A),  where 


A  =  Vu  +  (Vu)*  - 1  (V  •  u)I ,  (30) 

is  the  rate  of  deformation  tensor,  then  T  is  referred  to  as  a  viscous  stress  tensor.  Note,  that  for 
incompressible  flow,  the  divergence  of  the  velocity  is  zero.  However,  we  include  it  here  in  the  definition 
of  A  since  our  model  includes  compressible  (nonliquid)  regions.  For  a  Newtonian  fluid,  the  dependence 
of  the  stress  on  the  rate  of  deformation  tensor  is  linear,  (cf.  [24]) 

T  =  /t(4>)A,  (31) 

where  the  coefficient  fi  =  p(f)  may  depend  on  the  volume  fractions  of  the  fluids. 


Elastic-Plastic  Materials 


For  the  material  stress,  one  possibility  is  a  linear  elastic-plastic  model  in  the  form: 

f  C(E-H)  if\\E-U\\<EY 
|A(E-n)±f  iy |e - n||  =  £/ 

where  C  and  A  are  linear  operators  (tensors)  relating  the  strain  E  to  the  stress,  n  is  the  permanent  plastic 
deformation,  EY  is  the  yield  strain ,  f  is  the  yield  stress,  and  ||*||  is  an  appropriate  norm  acting  on  the 
strain  tensor.  The  strain  tensor  satisfies  the  system  of  equations  (this  represents  an  approximation  to  the 
strain) 


— +  u«VE  =  A, 

dt 


(33) 


and  n  is  evolved  using  the  equation 


dt 


+ u • vn  = 


<3E 

~dt 


+  u» 


VE 


,y||E-n||  =  Er 
,y||E-n|j<£/ 


(34) 


The  Eqs.  (33)  and  (34)  each  represent  an  additional  six  equations  (three  in  2-D  or  axi-symmetry)  for  the 
tensor  components  of  E  and  II.  Constitutive  laws  of  the  form  Eq.  (32)  depend  only  on  the  strain  of  the 
material  and  are  often  referred  to  as  “rate-independent.” 

Viscous  Models 

Alternatively,  it  may  be  more  appropriate  to  use  a  rate-dependent  material  description.  A  nonlinear 
rate  dependence  can  be  incorporated  directly  (without  the  need  for  introducing  extra  variables  as  with  Eq. 
(33))  into  the  momentum  equations  using  a  non-Newtonian  fluid  approximation.  One  non-Newtonian 
fluid  that  includes  plastic  behavior  is  the  Bingham  Model  [24,  25,  27].  This  model  uses  a  discontinuous 
stress-strain  rate  law  and  may  be  described  using  (cf.  Ref.  24,  p.  103) 
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T  = 


f 

V*(A:A) 


i/|(T:T)>f2, 


(35a) 


A  =0,  ifU  T:T)<f2. 


(35b) 


According  to  Bird  et  al.  (Ref.  25,  p.  227)  “The  viscoplastic  models  are  particularly  useful  for  describing 
liquids  with  large  amounts  of  suspended  solids,”  and  may  be  appropriate  for  the  limiting  case  of  saturated 
sand.  An  approximation  to  the  discontinuous  stress  described  in  Eq.  (35)  has  been  recently  proposed  by 
Papanastasiou  and  Boudouvis  in  Ref.  28.  This  exponential  model  has  the  form 


T=  A>  +  - 


l  -  exp|-a^  (A :  A) 


\'2 


HA:  A) 


A. 


(36) 


A  comparison  of  using  Eq.  (36)  with  various  a  compared  to  Eq.  (35)  is  shown  in  Fig.  1.  For  a  =  0  this 
model  reverts  to  Newtonian  viscosity.  It  approaches  the  Bingham  model  Eq.  (35)  as  a  — » 


Fig.  1  —  Stress-strain  rate  graphs  using  the  Papanastasiou  model  Eq.  (36) 


Viscoplastic  Flow  in  a  Circular  Tube 

Consider  the  steady-state  flow  of  a  viscous  fluid,  through  a  circular  tube  of  radius  R  and  length  L , 
with  pressures  P0  and  PL  specified  at  the  ends  of  the  tube.  In  this  case,  the  radial  velocity  component  is 
zero  and  the  axial  velocity  v  =  v(r)  depends  only  on  the  radial  location  r.  Both  the  deformation  and  stress 
tensors  have  a  single  non-zero  component.  By  considering  the  conservation  of  momentum  Eqs.  (28)  to 
(29)  including  the  hydrostatic  pressure  due  to  gravity  into  PL  and  integrating  in  space  along  both  the  axial 
and  radial  component,  yields 
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dv 

When  PL>P0,  it  follows  that  v'(r)  =  —  >0.  After  applying  Eq.  (36),  we  obtain 

dr 

T  =  jU0v'  +  f (1  -  ) .  (38) 

Using  a  no-slip  boundary  condition  v(/?)  =  0,  the  solution  to  Eq.  (37)  to  (38)  may  be  expressed  as 


j\ 

v(r)=-J 


g 


J  Pi-Po 


2  L 


ids , 


(39) 


^  _i  y 

where  g  '  denotes  the  inverse  of  g  with  g(y)  =  fl0y  +  f(l  -  e  "' ) .  In  the  case,  when  a  =  0,  g  (y)  =  —  and 

Mo 

the  usual  Newtonian  viscous  parabolic  flow  follows  directly  from  Eq.  (39): 


P  -P 
v(r)  =  —  \  L  0 


r 


l- 


4m/  j 

For  the  Bingham  model,  the  velocity  is  given  by  (see  Ref.  24) 


(r\ 

\RJ 


(40) 


v(r)  = 


rn-ni 

R2 

(  r  \ 
1-^ 

4 M/  j 

l  R) 

PlzI* 
V  4Mo L 


R‘ 


V 


R 


1  +  — ■ 


r  2  rn 


R  R 


r<rn 


r>r« 


(41) 


where  rn  = - . 

P  —P 

rL  ro 

Figure  2  shows  solutions  for  v  corresponding  to  different  values  of  the  parameter  a.  The  solutions  for 
0<a<°°  were  computed  using  Newton’s  method  to  determine  g“\  with  trapezoid  rule  approximation  of 
the  integral  Eq.  (39).  The  maximum  relative  error  in  velocity  between  the  Bingham  model  and  the 
exponential  model  Eq.  (36)  is  approximately  2 la  for  a  >10.  Note  the  “plug  flow”  behavior  of  the 
Bingham  model  for  r  <  r0. 

GENERALIZED  HYDRODYNAMICS  ALGORITHM  FOR  NON-NEWTONIAN  FLUIDS 

Here  we  describe  a  temporal  discretization  for  the  momentum  Eq.  (28)  in  the  case  of  a  non- 
Newtonian  fluid.  That  is,  we  will  consider  Eqs.  (28)  and  (29)  with  T  defined  in  Eq.  (36).  For  our  multiple 
species  case,  the  material  coefficients  fj.0  and  f  will  depend  on  the  volume  fractions.  To  clarify  the 
notation  for  the  functional  dependence  of  T  on  u  and  <|),  we  express 


T(<|>,  u)  =  /(4>,  A(u))A(u) , 


(42) 
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0  ro  R 


Fig.  2  —  Velocity  profiles  for  model  flows  in  a  circular  tube 


where  A  is  a  linear  operator  on  u  defined  in  Eq.  (30)  and 


/(4>.  A)  =  /l0  (<)>)  + 


T(<t>) 


1  -  exp^-a^(A:  A)j 


/f(A:A) 


(43a) 


An  exponential  dependence  was  used  for  the  viscous  and  yield  stress  terms  for  the  saturated  sand,  based 
on  an  empirical  data  fit  for  slurry  flows  [29].  For  example,  the  dependence  of  the  yield  stress  on  the 
volume  fraction  of  saturated  sand  is  given  by 


f  (<{))  =  f  (0S )  =  f 


exp(K^)-\ 

exp(Kz)-l 


(43b) 


where  kx  is  an  empirically  determined  constant,  and  (with  a  slight  misuse  of  notation)  associates  the 
constant  T  appearing  in  the  right-hand  side  of  Eq.  (43b)  with  f  (1) .  The  same  functional  form  is  used  for 
MoW- 


The  first  step  of  the  algorithm  is  the  explicit  convection  step.  This  represents  an  approximation  to  the 
solution  without  including  any  of  the  terms  on  the  right-hand  side  of  the  momentum  Eq.  (28).  That  is 


pu  =  (pu)n  -  tV  •  (puu)"+1/2 , 

^ =c -tv.  ton  r,/2, 

for  i  =  1, ... ,  m ,  and 

m 

P  =  IPto- 

i=l 


(44a) 

(44b) 

(44c) 


In  Eq.  (44)  the  superscript  rc  +  1/2  indicates  an  approximation  at  the  half-time  step.  This  step  is 
implemented  using  a  predictor-corrector  Godunov  scheme  (see  Refs.  11  and  13  for  details).  While  the 
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corrector  step  of  this  scheme  Eq.  (44)  ignores  the  pressure,  gravity  and  stresses  in  Eq.  (28),  these  terms 
are  included  in  the  predictor  step. 

The  redistribution  step  can  now  proceed  exactly  as  before  using  Eqs.  (10)  and  (11),  with  H  defined 
using  the  successive  constraint  approach  Eqs.  (24)  to  (27).  As  before,  the  redistributed  velocity  is  denoted 
u ,  and  the  redistributed  density  and  volume  fractions  are  simply  pw+1,  and  0"+1,  i  =  1, ... ,  m ,  respectively. 

The  viscous  stresses  are  included  in  the  diffusion  step  as  the  solution  u*  to 

p"+1  u*  -  pV  =  -rV  •  (puu)"+1/2  +  V2(Hu)  -  rpn+1/2gk  -  WPn~U2 

+tV  •  ((1  -  fl)T(<}>"+1/2,  u" )  +  0T(<t>"+1/2,  u* )) ,  (45) 

fj+l  _| _  Q11 

where  p 71+1/2  =  — —  —  (similarly  for  (j>),  and  6 is  a  temporal  discretization  parameter  for  the  diffusion 

term  0  <  0  <  1 .  The  value  0  =  0  corresponds  to  the  explicit  “Euler”  method,  0  =  1/2  to  the  second-order 
trapezoid  rule,  and  0  =  1  to  the  fully  implicit  “Backward  Euler”  method.  Note  that  except  for  the  lagging 
of  the  pressure  gradient  and  the  inclusion  of  the  momentum  redistribution  term,  this  equation  is  a  second- 
order  temporal  (trapezoid  rule)  discretization  of  the  momentum  Eq.  (28)  in  the  case  when  0  =  1/2.  Using 
Eqs.  (44)  and  (11)  with  Eq.  (45),  we  obtain 

p”+V  -  p"+1u  =  T0V  •  T((j)"+1/2, u*) -  r<9V  •  T((|)"+1/2,u"),  (46) 

where 

p"+1u  =  p"+1u  -  rpn+1,2gk  - 1 vr-1/2  +  tV  •  T(<t>”+1/2,  u").  (47) 

Subtracting  p"+1u”  and  setting  8  =  u*  -  u",  it  follows  from  Eqs.  (42)  and  (46)  that 

pn+1 5  -  T0V  •  (fWn+U2 ,  A(8  +  un ))  A(8))  =  pn+1  (u  -  un )  (48) 

+T0V  •  ({/(<T1/2,  A(8  +  u*))-  mn+l'2,  A(u"))}A(u”)). 

This  represents  a  nonlinear  equation  for  8.  Since  the  derivative  of /becomes  singular  as  a— and  | A| — >0, 
instead  of  solving  Eq.  (48)  using  Newton’s  method,  the  iteration  scheme 

p"+18m  -tOV •(f(f+U2 ,A(5lk-1]  +un))A(blk]))  =  pn+1(u-un)  (49) 

+T0V  •  ({/((T1/2,  A(8[*~1]  +  u"  ))-f(r+1,\A(un))}A(un )) 

is  used  with  8[0]  =  u  -  u”  for  k  >  1.  Each  step  of  Eq.  (49)  corresponds  to  the  solution  of  a  linear  system  of 
equations.  Upon  discretization  in  space,  the  linear  matrix  corresponding  to  the  left-hand  side  of  Eq.  (49)  is 
symmetric  and  positive  definite.  If /is  independent  of  A,  the  last  term  on  the  right-hand  side  of  Eq.  (49) 
vanishes  and  the  iterations  converge  in  a  single  step  (the  equation  is  linear  in  this  case).  In  the  special  case 
when  /  =  0,  Eq.  (49)  represents  a  trivial  identity  with  k  =  0.  Due  to  the  “stiff’  nature  of  the  linear 
equations  when  TOC  is  large,  the  value  0  =  1  was  used  for  all  the  computations  presented  here. 
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The  final  step  is  the  pressure  projection.  To  keep  this  part  of  the  algorithm  analogous  to  the  inviscid 
version,  the  lagged  pressure  gradient  must  be  added  back  to  the  momentum.  That  is,  the  projected 
momentum  is 


Pn+1u"+1  =  p"+v  +  r VP"-1'2  -  zVP"+1/2, 


(50) 


where  pn+u2  solves  the  equation 


with 


tV< 


_n+l 


VP' 


,n+l/2 


=  V  •  u*  in  D' 


,n+l 


(51) 


pn+,u*  =  p"+,u*  +tVP"'1/2.  (52) 

In  practice,  Eq.  (51)  is  not  solved  in  “delta”  form  since  the  domain  D"+1  and  the  Dirichlet  values  of  P  on 
the  boundaries  of  the  bubbles  change  with  time. 

NUMERICAL  RESULTS 

Three  distinct  numerical  problems  are  considered  in  this  section.  The  first  is  the  classical  problem  of 
viscous  flow  in  a  circular  tube  that  was  described  in  the  preceding  section  with  exact  solutions  plotted  in 
Fig.  2.  The  second  example  simulates  a  cylindrical  “slug”  of  fluid,  initially  at  rest  in  a  larger  cylindrical 
tank.  This  example  demonstrates  the  qualitative  capability  of  the  continuous  non-Newtonian  model  to 
simulate  viscoplastic  material  strength  in  the  fluid,  thereby,  allowing  “sand  castles”  to  stand.  Finally,  we 
compare  experiments  to  simulations  of  explosions  in  saturated  sand  under  a  layer  of  water. 


Computations  of  Flow  in  a  Circular  Tube 


For  our  first  example,  we  consider  the  problem  of  viscous  flow  in  a  circular  tube  (see  Fig.  2).  In 


particular,  we  have  considered  the  specific  problem  with  R  =  4 , 


Pl-Pq  , 
2  L  ’ 


jlIq  =  1 ,  and  f  =  2 ,  so  that  the 


plug-flow  region  for  the  limiting  Bingham  fluid  model  is  r  <  r0  =  2 .  The  convergence  of  the  computed 
solutions  to  the  “exact”  solutions  (using  Eq.  (38))  was  studied  by  decreasing  the  grid  spacing  in  the  radial 
direction.  In  particular,  a  sequence  of  grids  with  Nr  -  4,  8,  16,  and  32  uniform  cells  were  used.  In  each 
case,  only  four  uniform  cells  were  used  to  discretize  the  tube  along  its  length.  Initially,  the  velocity  inside 
the  tube  was  set  to  zero,  and  the  solutions  were  marched  in  time  to  steady  state,  driven  by  the  pressure 
boundary  conditions  at  the  tube  ends. 


Figure  3  shows  a  comparison  for  the  Newtonian  (a  =  0)  flow  case.  Even  on  the  coarsest  grid  with 
only  four  cells,  the  maximum  nodal  error  is  only  0.1%.  In  this  example,  the  solutions  for  all  the 
computations  are  close  to  the  exact  solution. 


Figure  4  shows  the  results  for  the  case  a  =  100,  which  is  within  1%  of  the  solution  of  the  limiting 
Bingham  plastic  case.  In  this  case,  the  computational  errors  were  much  larger  than  for  the  Newtonian 
case.  However,  it  is  apparent  from  Fig.  4  that  the  solutions  improve  as  the  grid  is  refined. 
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Fig.  3  —  Computed  and  exact  solution  to  Newtonian  flow  in  a  circular  tube 


Fig.  4  —  Computed  and  exact  solutions  to  generalized  Bingham  flow  (a  =  100)  in  a  circular  tube 


Table  1  lists  a  summary  of  the  computed  maximum  errors  and  rates  of  convergence.  As  noted  in 
Figs.  3  and  4  above,  the  errors  are  much  larger  for  the  a  =  100  case.  For  the  Newtonian  model,  the  nodal 
rate  of  convergence  was  about  3,  indicative  of  super-convergence  at  the  nodal  points.  In  general,  only  a 
second-order  convergence  rate  could  be  expected.  However,  for  the  Bingham  model,  the  rate  of 
convergence  was  only  about  1.5,  possibly  due  to  the  discontinuity  in  the  second  derivative  of  the  solution 
at  the  point  where  the  flow  becomes  plastic. 

Slug  Flow 

In  our  next  example,  we  examine  the  qualitative  behavior  of  a  model  of  a  Bingham  plastic  slurry.  It  is 
obvious  that  sand  (particularly  wet  sand)  does  not  flow  like  a  Newtonian  fluid  and  can  sustain  nontrivial 
deformations.  We  demonstrate  numerically,  how  sand  castles  can  be  simulated  using  the  Papanastasiou 
approximation  Eq.  (36)  to  the  Bingham  plastic  model. 
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Table  1  —  Errors  and  Rates  of  Convergence  for  Circular  Tube  Flow  Problem 


a  =  0 

Q 

11 

o 

o 

Nr 

K 

Rel. 

Rate 

K 

Rel. 

Rate 

Error% 

Error% 

4 

1.0 

1.06e-l 

0.4 

36.26 

8 

1.0 

1.27e-2 

3.06 

0.4 

13.72 

1.40 

16 

1.0 

1.63e-3 

2.96 

0.4 

5.09 

1.43 

32 

1.0 

2.50e-4 

2.70 

0.4 

1.76 

1.53 

Consider  a  cylindrical  slug  of  a  slurry  initially  at  rest  in  a  larger  cylindrical  tank  with  rigid  walls  in 
which  gravity  is  acting.  The  problem  has  been  nondimensionalized  by  setting  ps  =  1,  g  =  1,  PA  =  1,  and 
p  —  l  inside  a  cylinder  of  radius  1/4  and  height  1.  The  boundary  walls  of  the  cylindrical  tank  have  radius 
1/2.  A  uniform  grid  of  square  cells  of  size  h  —  \/37.  was  used  for  the  simulations.  For  comparison,  we 
repeated  the  simulation  for  several  different  cases  listed  in  Table  2. 


Table  2  —  Bingham  Parameters  for 
Different  Slug  Flow  Problems 


Case 

mm 

f 

a 

T 

A 

0 

0 

5.0e-4 

B 

0 

5.0e-4 

C 

0 

1 

50 

5.0e-4 

D 

0 

Km 

200 

5.0e-4 

E 

0 

200 

5.0e-4 

F 

0 

200 

5.0e-4 

G 

0 

200 

2.5e-4 

Case  A  is  the  flow  of  an  incompressible  inviscid  liquid  cylinder  and  is  included  for  the  sake  of 
comparison.  Figure  5  displays  density  profiles  and  velocity  vectors  for  this  case  at  various  times.  The 
density  contours  are  drawn  at  five  levels  p  =  <ps  =  jO.  1,0.3, 0.5, 0.7, 0.9}.  To  improve  visualization,  velocity 

vectors  are  drawn  only  from  the  centers  of  the  lower  left  comer  cell  within  each  2  by  2  block  of  cells.  A 
vector  is  drawn  whenever  p  >  0.01  to  highlight  spray  regions.  A  relative  scale  for  the  length  of  the 
velocity  vectors  is  used.  The  value  listed  below  each  frame  next  to  the  arrow  is  the  maximum  speed  in  the 
region  where  p>0.01  and  corresponds  to  the  length  of  the  adjacent  arrow.  The  second  value  in 
parentheses  is  the  maximum  speed  in  the  “liquid”  region  where  p  >0.9.  Note  the  spray  region,  which 
forms  as  the  liquid  falls  after  rising  up  the  side  walls.  The  sloshing  of  the  liquid  continues  for  a  substantial 
time  until  numerical  dissipation  finally  causes  the  liquid  to  come  to  a  level  rest  (not  shown). 

Case  B  represents  slug  flow  for  a  Newtonian  viscous  fluid  and  is  displayed  in  Fig.  6.  In  this  case,  the 
rising  of  the  liquid  up  the  side  walls  is  not  nearly  as  high  as  the  inviscid  case,  and  no  spray  is  formed.  The 
convergence  of  this  flow  to  the  trivial  steady  state  level  solution  is  evidenced  in  the  last  three  frames. 

Cases  C  and  D  are  approximations  to  non-Newtonian  Bingham  flow.  The  value  f  =  0.2  exceeds  the 
stress  due  to  gravity  at  the  top  of  the  slug  but  not  at  the  bottom.  Thus,  the  cylindrical  slug  will  flow 
similar  to  a  stagnation  flow  at  the  bottom,  while  the  top  region  will  fall  as  a  rigid  plug.  When  the  fluid  at 
the  bottom  reaches  the  cylindrical  walls,  the  stress  is  reduced  below  the  yield  stress  and  the  fluid  stops. 
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Fig.  5  —  Slug  flow  for  an  inviscid  incompressible  liquid  (Case  A) 


Figure  7  shows  the  results  for  the  Papanastasiou  approximation  to  Bingham  flow,  with  the  parameter 
a  =  50,  and  Fig.  8  displays  the  approximation  with  a  =  200.  Note  that  the  higher  value  for  a  causes  the 
fluid  in  the  region  where  the  yield  stress  is  exceeded  to  flow  at  an  earlier  time.  This  is  due  to  the  fact  that 
if  the  yield  stress  is  exceeded,  the  “effective  viscosity”  (the  slope  of  the  curves  depicted  in  Fig.  1) 
converges  to  fi0(=  0)  as  a  ->  Therefore,  the  region  of  yielding  flow  in  Case  C  experiences  a  higher 
viscosity  than  in  Case  D.  Note,  however,  that  the  final  steady-state  solutions  (T>4.00)  for  both  cases 
have  very  similar  nontrivial  profiles.  Next,  we  consider  the  case  when  the  spread  of  the  slug  does  not 
reach  the  wall  boundary.  Figure  9  displays  this  case  and  clearly  demonstrates  the  model’s  ability  to 
simulate  sand  castles. 

Figure  10  displays  a  slug  which  behaves  as  a  rigid  body.  The  slight  motion  as  indicated  by  a  small 
indentation  in  the  top  of  the  slug  and  a  small  contour  line  near  the  bottom  center  are  due  to  numerical 
error.  Fig.  11  shows  the  same  example  computed  with  the  time  steps  halved.  This  example  indicates  that 
while  small  time  steps  are  not  required  for  the  overall  stability  of  the  algorithm,  they  are  still  needed  for 
accuracy  when  the  Bingham  yield  stress  is  high. 

Explosive  Channeling  Results 

The  motivation  for  this  research  was  the  ability  to  predict  channels  and  craters  formed  in  sand  from 
underwater  explosions.  A  small-scale  2-D  test  box  has  been  constructed  at  the  University  of  Maryland  to 
visually  examine  the  dynamics  of  explosive  channeling  [33].  This  test  box  (referred  to  as  the  UMD  2-D 
Box)  is  depicted  with  dimensions  in  Fig.  12.  In  this  box,  a  layer  of  sand  and  water  are  sandwiched 
between  a  plexiglass  front  and  aluminum  plate  back.  The  distance  between  the  front  and  back  plates  is 
small  with  respect  to  the  size  of  a  free-field  spherical  bubble  in  water  formed  from  the  detonation. 
Therefore,  neglecting  viscous  effects  at  the  plate  boundaries,  the  dynamics  of  the  bubble,  water,  and  sand 
will  be  primarily  2-D. 
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Fig.  6  —  Slug  flow  for  an  viscous  incompressible  liquid  (Case  B) 
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Fig.  7  —  Slug  flow  for  Case  C  (/Iq  =  0,t  =  0.2,  and  a  =  50) 


T=0.00 


T- 1 .00 


0.00  (0.00)  -0.01  (0.01) 


T-2.00 


T  =  3.QQ 


T  =  4.00 


fry  ‘ 

-1.16  (0.37) 


Fig.  8  —  Slug  flow  for  Case  D  (jU0  =  0,t  =  0.2,  and  a  =  200) 
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Fig.  9  —  Slug  flow  for  Case  E  (p0  =  0,T  =  0.3,  and  a  =  200) 


Fig.  1 0  —  Slug  flow  for  Case  F  (n0  =  0,  i  =  0.5, a  =  200,  and  r  =  5e  -  4) 


Fig.  1 1  —  Slug  flow  for  Case  G  (/i0  =  0,  t  =  0.5,  a  =  200,  and  t  =  2.5e  -  4) 


The  charges  used  in  the  experiments  were  RP-80  detonators  containing  0.203  grams  of  explosive. 
The  experiments  were  photographed  using  a  digital  camera  which  recorded  500  frames  per  second  at  a 
resolution  of  512  x  512  pixels  and  a  depth  of  8-bits  (grayscale).  (Additional  details  of  the  experimental 
setup  and  instrumentation  may  be  found  in  Ref.  33.)  The  images  shown  were  enhanced  by  first  reducing 
noise  with  a  Total  Variation  Minimization  algorithm  [30],  followed  by  contrast  equalization  [31].  The 
noise  reduction  was  sufficiently  small  (the  variance  of  the  image  was  reduced  by  only  2%),  so  that  no 
important  features  were  removed. 

Figure  13  contains  images  from  an  experiment  in  which  the  initial  height  of  water  Hw  was  2.5  in., 
and  the  depth  of  burial  Db  was  approximately  1  in.  (measured  from  the  top  of  the  sand  to  the  center  of  the 
charge).  In  these  images,  the  sand  was  “striped”  using  dye  to  improve  the  visualization  of  the  channel 
formation.  Each  stripe  is  approximately  1  in.  thick.  Shortly  after  detonation,  a  gas  bubble  forms  and 
continues  to  expand  until  it  attains  it  maximum  size  at  t  —  0.008.  During  the  expansion,  a  water  column 
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Fig.  12  —  The  UMD  2-D  box  experimental  setup 

(unfortunately  mostly  obscured  by  the  ruler)  rises  above  the  bubble.  The  outline  of  the  bubble  is  partially 
obscured  by  the  light  sand,  particularly  on  the  sides  above  the  initial  sand  level.  The  bubble  is  best 
observed  at  t  =  0.010,  when  it  is  marked  by  dark  spots  surrounded  by  a  mixture  of  light  sand  and  water. 
After  the  bubble  maximum  is  attained,  low  pressure  inside  the  bubble  causes  it  to  contract  as  seen  in  the 
images  corresponding  to  times  0.010  <t  <0.016.  The  important  feature  to  note  in  these  images  is  the 
formation  of  a  downward  water  jet  moving  through  the  center  of  the  bubble.  This  jet  appears  as  a  dark 
object  near  the  top  of  the  bubble  at  t  =  0.012,  which  moves  downward  at  t  =  0.014  and  t  =  0.016.  The 
channel  walls  appear  to  be  widened  after  the  jet  impacts  the  bottom  and  flows  up  the  sides  during 
0.018  <  t  <  0.100.  The  sand  is  suspended  in  the  slurry  mixture  at  t  =  0.200,  and  t  =  0.500  and  begins  to 
settle  out  and  flow  back  into  the  channel  at  t  =  1.000.  After  t  =  1.500,  the  channel  reaches  a  steady  state, 
where  the  true  crater  has  been  partially  refilled  by  sediment  to  form  a  shallower  apparent  crater.  The 
observation  of  “true”  and  “apparent”  craters  also  appears  from  experiments  of  explosives  buried  in  soil 
without  a  water  layer,  as  discussed  in  Ref.  1.  The  use  of  dyed  sand  in  the  experiments  enhances  the 
distinction  between  the  true  and  apparent  craters  at  t  =  1.500. 

Computations  were  performed  under  the  assumption  that  the  experiment  is  nearly  2-D.  Empirical 
laws  are  used  to  initialize  the  computation  based  on  the  theory  described  in  Refs.  5  and  6.  For  a  2-D 
computation,  the  explosion  is  initialized  as  a  high-pressure  cylinder  of  gas  in  fluids  which  are  initially  at 
rest.  Effects  of  the  shock  wave  are  neglected  in  this  study.  The  bubble  pressure  is  assumed  to  behave 
adiabatically  PB  =  CVBy ;  where  C  is  a  constant,  VB  is  the  bubble  volume  and  y  is  the  ratio  of  specific 
heats  of  the  bubble  gases.  Following  Refs.  5  and  6  and  using  the  superscript  (2D)  to  denote  a  2-D  value 
(no  superscript  is  used  for  3-D),  we  have 


at  = 

(53) 

“a°'=b'n' 

(54) 

a(2D)  =  aV2, 

(55) 

r>4^3/!. 

(56) 
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(57) 

(58) 


where  A™  is  the  initial  (minimum)  radius  of  the  cylinder,  qaD)  is  the  empirical  minimum  radius 
constant  depending  on  the  charge  type,  and  M  is  the  mass  per  unit  length  of  the  charge.  The  free-field 
maximum  radius  is  A™  and  a(2D)  =  A™  /  A™-  The  other  empirical  constant  depending  on  the  charge 
type  is  /i2£>),  often  referred  to  as  the  maximum  radius  constant.  As  was  discussed  in  detail  in  Ref.  5,  Eq. 
(57)  is  often  approximated  using 


(2D)  _  r(2D) 


M 


111 


Pj 


12 


(59) 


which  is  a  valid  only  when  the  charge  is  very  deep  relative  to  the  maximum  bubble  radius  Instead  we 
use  the  corrected  Eq.  (57),  which  requires  the  solution  of  a  scalar  nonlinear  equation  for  the  determination 
of  ai2D)  and,  hence,  PB°.  The  traditional  units  used  for  underwater  explosion  bubbles  are  feet  for  length, 
pounds  for  charge  weight,  and  feet  of  water  for  the  pressure.  The  pressure  at  infinity  P  is  taken  to  be  the 
hydrostatic  head  at  the  initial  charge  depth.  For  example,  under  ambient  conditions  in  fresh  water,  we 
have 

=33.9  +  d,  (6°) 


where  d  is  the  initial  depth  in  feet,  and  P„  is  measured  in  units  of  feet  of  water.  Initial  conditionsjor  t  e 
computations  shown  here  were  derived  from  the  above  equations  with  constants  J„  =  11.0,  q  =  0.272 
and  7  =  1.3.  This  value  used  for  J _  is  approximately  30%  lower  than  expected  for  RP-80  detonators.  It  is 
conjectured  that  the  smaller  than  expected  bubble  sizes  are  due  to  either  the  lack  of  rigidity  of  the  front 
and  back  faces  of  the  2-D  box  or  the  rigid  body  motion  of  the  entire  2-D  box  during  the  tests.  The  entire 
box  would  typically  “jump”  about  1  in.  after  the  charge  was  detonated. 

Since  the  radius  ratio  a(2D)  =  41.6  is  large,  and  it  is  desirable  to  resolve  the  bubble  at  both  its 
minimum  and  maximum  volumes;  different  grids  were  used  at  different  time  intervals  of  the  simulation. 
All  grids  contained  a  rectangular  region  0<*<X*  and  ZB  <z<Zr  of  square  cells  of  size  h.  This 
rectangle  is  surrounded  by  coarser  cells  that  are  stretched  to  the  physical  boundaries  of  the  2-D  box  _ We 
use  the  convention  that  *  =  0  is  a  symmetry  line  and  z  =  0  is  the  initial  location  of  the  air-water  interface. 
The  total  number  of  cells  in  the  entire  domain  is  denoted  by  NxxNz.  Table  3  shows  a  description  of  the 
grids  and  time  intervals  in  which  they  were  used.  Grid  1  is  used  to  resolve  the  very  small  bubble  at  early 
times  until  it  grows  to  approximately  five  times  its  initial  radius.  The  flow  field  is  then  mapped  onto 
Grid  2,  such  that  mass  and  momentum  are  conserved  and  the  computation  is  restarted.  Grid  3  is  a  coarse 
grid  used  for  the  long-time  wave  dynamics  after  the  bubble  has  vented. 
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Table  3  —  Grid  Parameters  Used  for  the  UMD  2-D  Box  Simulations 


Grid 

HI 

■SI 

1 

0.0  <t  <0.00025 

0.0012 

0.6 

-4.3 

-3.1 

70x150 

2 

0.00025  <  t  <  0.2 

0.012 

9.0 

-8.4 

8.4 

76x184 

3 

0.2  <t  <1.5 

0.024 

9.0 

-8.4 

8.4 

38x97 

Figure  14  shows  a  simulation  of  the  UMD  2-D  box  experiment  using  the  model  described  in  this 
report.  For  this  computation,  the  yield  stress  (nondimensionalized  by  dividing  by  the  ambient  pressure) 
was  f  =  0.59.  For  these  simulations,  the  ratio  of  yield  stress  to  plastic  viscosity  was  f //U0  =83.  The 
parameter  ct  =  10  was  used  in  the  Papanastasiou  model  Eq.  (36).  The  early  expansion  of  the  bubble,  jet 
formation  and  collapse  appear  to  be  reasonably  well  predicted.  However,  the  effect  of  the  jet  impact  on 
the  bottom  of  the  channel  is  much  more  severe  in  the  model  than  in  the  experiment.  Furthermore,  the 
model  is  not  able  to  predict  the  washback  or  sedimentation  of  the  sand  back  into  the  channel.  The  final 
(steady-state)  channel  (t  >1.0)  is  substantially  deeper  than  either  the  apparent  or  true  craters  formed  from 
the  experiment.  Also,  the  “lip”  (area  of  raised  sand  at  the  channel  edges)  is  higher  and  narrower  than  in 
the  experiment. 

Figure  15  shows  a  simulation  of  the  same  experiment  but  with  a  scaled  yield  stress  of  f  =  0.059.  In 
this  case,  the  steady-state  channel  is  in  reasonably  good  agreement  with  the  experimental  channel. 
Unfortunately,  when  comparing  Fig.  13  to  Fig  15,  one  sees  little  resemblance  of  the  dynamics  from  the 
simulation  to  the  experiment.  With  the  lower  yield  stress,  the  bubble  jets  at  a  much  later  time  than  in  the 
experiment.  Also,  the  simulation  shows  the  bubble  migrating  below  the  initial  charge  location.  There  is 
no  evidence  of  this  in  the  images  from  the  experiment.  This  example  clearly  shows  the  danger  in  simply 
comparing  the  final  results  of  a  simulation  to  an  experiment,  particularly  when  “violent”  nonlinear 
hydrodynamics  is  occurring. 

Figure  16  compares  the  final  channels  from  the  experiment  and  simulations  using  various  yield 
stresses.  It  is  evident  that  none  of  the  simulations  are  capable  of  providing  an  accurate  description  of  the 
final  channel  profile.  With  increasing  yield  stresses,  the  depth  of  the  channel  approaches  that  of  the 
experiment.  However,  the  width  of  the  channel  also  decreases  and  becomes  less  like  the  experiment. 

Based  on  Fig.  16,  the  current  model  does  not  appear  to  be  capable  of  reliably  predicting  the  complex 
interaction  between  the  sand  and  water  upon  the  bubble  jet  collapse  or  the  later  sedimentation  of  the  sand 
back  into  the  channel.  However,  the  early  time  bubble  formation  can  be  accurately  predicted.  In  another 
experiment  in  the  UMD  2-D  box,  in  which  Hw  =  3.0  in.  and  Db  =0.6  in.,  the  outline  of  the  bubble  at  its 
maximum  volume  was  clearly  displayed  in  the  video  image.  This  image  was  compared  to  simulations 
using  various  yield  stresses  in  Fig.  17.  In  this  figure,  the  agreement  of  the  simulation  with  f  =  0.295  to 
the  experiment  is  remarkable.  In  the  case  of  no  yield  stress  (here  the  saturated  sand  is  treated  as  a  dense 
incompressible  liquid),  the  bubble  has  an  oval  shape  and  is  much  deeper  in  the  sand  than  the  experiment 
indicates.  Incorporating  a  yield  stress  produces  more  realistic  bubble  profiles  in  which  the  bubble  radius  is 
larger  in  the  water  region  than  the  sand.  This  is  intuitively  correct  and  verified  by  the  image  from  the 
experiment  in  Fig.  17. 
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Pig  17  —  Comparison  of  bubble  at  maximum  volume  for  the  case  Hw  -  3.0  in.  and  Db  -  0.6  in. 


SUMMARY  AND  CONCLUSIONS 

In  an  attempt  to  model  explosive  channeling  in  water-covered  sand,  we  have  generalized  a  model 
which  has  been  successful  for  violent  free-surface  motions  with  a  single  incompressible  liquid  species. 
The  new  model  included  viscoplastic  forces  for  saturated  sand,  so  that  nontrivial  steady-state  conditions 
(e.g.,  channels)  could  be  predicted.  Due  to  the  stiff  nature  of  the  equations  employed  to  approximate  a 
Bingham  plastic,  a  fully  implicit  numerical  treatment  was  required  for  stability.  The  model  was 
successfully  tested  on  a  benchmark  problem  of  slug  flow  in  a  circular  cylinder,  for  which  the  exact 
solution  is  known.  It  was  also  tested  on  a  problem  of  a  free-standing  cylinder  of  material  under  the 
influence  of  gravity,  which  qualitatively  demonstrated  the  capability  to  predict  nontrivial  steady-state 
deformations. 

For  simulations  of  channeling,  the  early  time  bubble  expansion  and  first  collapse  were  accurately 
predicted,  provided  that  an  appropriate  value  for  the  yield  stress  was  used.  This  result  is  encouraging  and 
suggests  that  the  Bingham  model  is  appropriate,  at  least  during  the  first  period  of  the  bubble  dynamics. 
The  effect  of  the  water  jet  impacting  the  bottom  of  the  channel  appeared  to  be  overpredicted  in  the 
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simulations.  Also,  no  mechanism  for  the  sedimentation  or  the  wash  back  of  the  sand  refilling  the  channel 
is  provided  in  the  model.  This  wash  back  is  responsible  for  the  substantial  difference  observed  between 
the  apparent  and  true  channels. 

So  far,  we  have  considered  the  treatment  of  multiple  species  with  only  one  flow  field  representing  the 
velocity  of  the  mixture  of  the  species.  This  formulation  does  not  model  the  flow  of  water  through  the 
sand,  which  can  influence  the  sand  strength  and  describe  the  settling  of  sand  mixed  with  water.  In  future 
work,  we  will  describe  a  model  that  includes  a  separate  velocity  field  for  each  species.  This  model  will 
also  account  for  the  porosity  of  sand,  by  constraining  its  volume  fraction  to  be  bounded  by  a  quantity  less 
than  unity.  The  permeability  of  sand  can  be  modeled  by  the  inclusion  of  drag  terms  in  the  momentum 
equations.  Finally,  more  physically  realistic  stress  models  for  the  saturated  sand,  such  as  those  satisfying 
the  Mohr-Columb  criterion  [  1 9-2 1  ]  or  Prager-Drucker  condition  [22,  23]  will  be  examined. 
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